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ABSTRACT 

A  variational  data  assimilation  system  is  developed  for  the  stationary,  homogeneous  portion  of  the  wave 
model  Simulating  Waves  Nearshore  (SWAN).  The  system  is  based  on  a  numerical  adjoint  constructed  for  the 
discrete  forward  SWAN  code;  its  performance  is  compared  to  that  of  an  earlier  system  based  on  a  discretized 
analytical  adjoint  (Walker;  Veeramony  et  al.).  This  paper  describes  the  development  and  validation  of  in¬ 
dividual  numerical  adjoint  subroutines,  followed  by  the  testing  and  evaluation  of  the  assimilation  system  as 
a  whole  with  an  idealized  twin  experiment  and  with  data  from  Duck,  North  Carolina.  In  the  twin  experiment, 
the  present  system  performs  on  par  with  that  of  Walker.  Estimates  of  wave  spectra  and  spectral  statistics  also 
compare  well  to  measured  spectral  data  at  Duck,  North  Carolina.  The  error  in  these  estimates  is  partly  due  to 
the  exclusion  of  nonlinear  source  and  sink  terms  from  the  adjoint  and  partly  due  to  different  spectral  pro¬ 
cessing  techniques  used  for  different  types  of  instruments. 


1.  Introduction 

The  wave  model  Simulating  Waves  Nearshore  (SWAN) 
(Booij  et  al.  1999)  solves  the  spectral  action  balance 
equation  [Eq.  (1)]  to  produce  nearshore  wave  forecasts  and 
climatologies.  It  is  widely  used  by  the  coastal  modeling 
community  and  is  part  of  a  variety  of  coupled  ocean-wave- 
atmosphere  model  systems.  In  forecasting  wave  conditions 
for  specific  locations  or  events,  boundary  conditions  for  the 
local  SWAN  domains  are  generally  obtained  from  regional 
or  global  simulations  with  WAVEWATCH  III  (Tolman 
2009)  and  the  Wave  Model  (WAM;  WAMDI  1988).  The 
accuracy  of  nearshore  wave  estimates  is  highly  dependent 
on  the  quality  of  these  boundary  inputs.  Even  small  errors 
in  wave  energy,  period,  and  direction  at  the  boundary  can 
grow  significantly  as  SWAN  propagates  a  wave  field  into 
shallow  water.  With  a  system  that  utilizes  measured  data  to 
correct  the  boundary  condition  for  the  SWAN  model,  such 
errors  can  be  considerably  reduced. 

The  principal  spectral  action  balance  equation  in 
SWAN  may  be  expressed  as 

dL  +  V.(CN)  =  S^,  (1) 

dt  a 
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where  N  is  the  spectral  action  density,  obtained  by 
dividing  the  spectral  energy  density  by  an  intrinsic 
representative  wave  frequency  a.  The  first  term  of 
Eq.  (1)  expresses  the  rate  at  which  the  action  density  is 
changing  with  time.  In  the  second  term,  the  propagation 
of  action  density  in  both  physical  (x,  y )  and  spectral  (cr,  6) 
space  is  represented  by  the  dot  product  of  vector  gra¬ 
dient  V  =  (d/dx,  d/dy,  d/dcr,  d/dd)  with  vector  velocity 
C  =  (Cx,Cy,Ca,Ce)  and  action  density  N.  The  vector 
C  represents  the  propagation  rate  of  wave  energy  and 
includes  contributions  from  both  wave  group  velocity 
and  ambient  current.  On  the  right-hand  side  of  Eq.  (1), 
St*  *s  a  consolidated  spectral  source  term  that  may 
include  contributions  from  nonlinear  wave-wave  in¬ 
teractions,  wind  wave  forcing,  dissipation  due  to 
breaking,  bottom  friction,  and  whitecapping  [see  Booij 
et  al.  (1999)  for  additional  details].  Once  initialized  with 
bathymetry,  boundary  wave  data,  initial  conditions,  and 
source  term  configuration,  SWAN  solves  Eq.  (1)  to  pre¬ 
dict  the  evolution  of  the  wave  spectrum  throughout  the 
modeled  coastal  region. 

While  data  assimilation  has  been  a  component  of  at¬ 
mospheric  modeling  since  the  1960s,  it  has  only  recently 
been  introduced  in  ocean  wave  models.  A  number  of 
nearshore  studies  have  used  assimilated  wave  data  from 
in  situ  instruments  or  video  to  estimate  a  variety  of 
nonwave  nearshore  parameters,  including  bottom  ba¬ 
thymetry  (Lippmann  and  Holman  1990;  van  Dongeren 
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et  al.  2008;  Narayanan  et  al.  2004;  Aarninkhof  et  al. 
2005),  alongshore  currents  (Chickadel  et  al.  2003),  and 
bottom  friction  (Keen  et  al.  2007).  Early  attempts  to 
improve  wave  predictions  with  data  assimilation  used 
relatively  simple  optimal  interpolation  methods  with 
observed  wave  heights  and  periods  to  nudge  model  re¬ 
sults  toward  more  accurate  values  (Bauer  et  al.  1992). 
More  recent  efforts  have  improved  model  performance 
further  by  incorporating  wave  spectral  parameters 
(Hasselmann  et  al.  1997;  Aouf  et  al.  2006;  Portilla  2009). 
However,  these  features  have  primarily  been  used  in 
larger-scale  models  such  as  WAVEWATCH  III  and 
WAM.  The  public  domain  version  of  SWAN  does  not 
yet  include  data  assimilation  in  any  form. 

Walker  (2006)  developed  a  data  assimilation  system 
that  used  measured  data  from  inside  the  computa¬ 
tional  domain  to  correct  offshore  boundary  inputs  for 
SWAN.  In  this  system,  the  adjoint  of  reduced  SWAN 
equations  was  derived  analytically  and  then  discretized. 
This  “analytical  adjoint”  was  used  to  propagate  the 
model-observation  errors  to  the  boundary,  where  only 
stationary  (i.e.,  steady  state)  conditions  were  considered, 
and  all  nonlinear  source  and  sink  contributions  were  set 
to  zero,  as  in 

V  •  (CAT)  =  0.  (2) 

This  approach  assumes  that  energy-generating  or  dissi¬ 
pating  effects  like  wind  forcing,  wave  breaking,  and 
bottom  friction  are  negligible.  The  steady  state,  linear 
form  of  SWAN  represented  by  Eq.  (2)  is  commonly  used 
to  investigate  longer-term  effects  of  an  average  or  rel¬ 
atively  constant  wave  climate.  Recent  validation  tests 
indicate  that  Walker’s  assimilation  system  leads  to  im¬ 
proved  predictions  of  integrated  spectral  wave  proper¬ 
ties  like  significant  wave  height,  peak  period,  and  mean 
direction  in  deeper  water,  as  well  as  limited  improve¬ 
ments  in  estimates  of  the  full  2D  spectrum  (Veeramony 
et  al.  2010).  However,  the  technique  does  not  capture 
changes  to  the  spectra  that  occur  due  to  the  contribu¬ 
tions  of  nonlinear  source  terms  described  above.  Since 
the  algorithm  only  corrects  boundary  conditions  off¬ 
shore,  an  ad  hoc  treatment  of  lateral  boundaries  is  re¬ 
quired.  These  shortcomings  generally  lead  to  significant 
errors  in  predicted  wave  heights  and  currents  in  shallow 
water. 

The  analytical  adjoint  approach  cannot  be  adapted  to 
incorporate  nonlinear  wave  transformations  in  a  consis¬ 
tent  manner.  It  has  been  demonstrated  (Biicker  et  al. 
2011;  Giresse  and  Walther  2004)  that  discretizing  the 
analytical  adjoint  model  (as  done  in  Walker  2006)  leads 
to  an  erroneous  numerical  gradient  of  the  cost  function. 
As  a  result,  an  analytical  adjoint  system  generally  will 


not  pass  the  asymptotic  gradient  test  of  its  cost  function 
(cf.  Jarvinen  1998).  The  proper  and  consistent  way  to 
obtain  the  numerical  gradient  of  the  cost  function  is  to 
derive  the  adjoint  of  the  discretized  forward  model 
(hereafter  called  the  “numerical  adjoint”).  The  differ¬ 
ence  in  the  numerical  gradients  obtained  from  the  two 
adjoint  types  is  small  when  working  with  a  stationary 
linear  system,  but  it  can  be  significantly  larger  for  adjoints 
to  nonstationary  and/or  nonlinear  models  (Favennec 
2005). 

For  environments  where  wind  forcing,  bottom  friction, 
currents,  or  triad-quadruplet  interactions  are  significant, 
effective  adjoint-based  assimilation  of  wave  data  thus 
requires  a  numerical  approach.  A  numerical  adjoint  to 
SWAN  can  include  adjoint  subroutines  for  all  nonlinear 
sources  and  sinks  of  wave  energy  in  the  model,  regard¬ 
less  of  their  initial  discrete  form.  An  assimilation  system 
featuring  a  complete  numerical  adjoint  will  enable  SWAN 
users  to  improve  spectral  estimates  with  assimilated  wave 
data  from  nearly  any  domain,  including  storm-generated 
chop  and  swell,  shallow  surf  zones,  and  other  nonlinear 
wave  environments.  Construction  of  such  an  adjoint  is  la¬ 
bor  intensive,  however,  requiring  each  adjoint  subroutine 
to  be  created  and  validated  separately. 

The  present  paper  describes  the  development  and 
testing  of  a  variational  data  assimilation  system  based 
upon  a  numerical  adjoint  that  is  limited  to  the  stationary, 
linear  SWAN-governing  equations  [Eq.  (2)].  The  as¬ 
similation  system  described  herein  is  suboptimal,  im¬ 
plementing  strong-constraint  variational  assimilation  in 
which  just  the  boundary  condition  is  controlled.  Although 
fully  consistent  four-dimensional  variational  data  assimi¬ 
lation  (4DVAR)  systems  have  been  built  and  validated  for 
large-scale  circulation  models,  such  as  the  Navy  Coastal 
Ocean  Model  (NCOM;  Ngodock  and  Carrier  2013),  the 
Ocean  Parallelise  (OP A)  model  (Weaver  et  al.  2003),  and 
the  Regional  Ocean  Model  System  (ROMS;  Moore  et  al. 
2004),  this  is  the  first  time  one  has  been  created  for  a  wave 
model.  The  purpose  of  this  study  is  to  introduce  this  new 
approach  to  the  nearshore  community,  demonstrate  its 
feasibility,  and  evaluate  the  initial  suboptimal  system 
through  a  direct  comparison  with  the  equivalently  re¬ 
stricted  analytical  adjoint  created  by  Walker  (2006) 
and  tested  by  Veeramony  et  al.  (2010).  Consequently, 
the  same  physical  assumptions  are  adopted  here  as  in 
Walker’s  analytical  approach  (i.e.,  operating  only  under 
stationary  conditions  and  excluding  nonlinear  sources 
and  sinks  in  the  adjoint).  The  following  two  sections 
summarize  the  general  structure  and  theory  of  the  as¬ 
similation  system,  the  development  of  the  adjoint  to  the 
simplified  model,  and  the  initial  validation  tests.  In  sec¬ 
tion  4,  the  new  system’s  performance  is  compared  directly 
to  that  of  the  Walker  (2006)  system  in  a  set  of  twin 
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Fig.  1.  Assimilation  system  flowchart. 


experiment  simulations,  and  then  further  examined 
using  measured  nearshore  spectral  data  and  bathy¬ 
metry  from  Duck,  North  Carolina.  Additional  discussion 
is  offered  in  section  5,  followed  by  a  brief  summary  of 
overall  conclusions. 

2.  Assimilation  system 

As  mentioned  above,  in  the  assimilation  system  pre¬ 
sented  here,  the  numerical  adjoint  encompasses  only  the 
stationary,  homogeneous  part  of  the  SWAN  model,  in¬ 
cluding  linear  subroutines  for  the  transport  and  re¬ 
fraction  of  wave  action  but  no  sources  or  sinks.  The  data 
assimilation  procedure  is  built  around  an  objective  cost 
function  J(N ),  which  is  used  to  measure  the  aggregate 
model  and  measurement  error  at  data  locations  within 
the  model  domain  and  has  the  following  form: 

J(N)  =  (BC  -  BC6)tQb£(BC  -  BC6) 

+  £  O'*  -  "A)TR  '(>>  -  HkN).  (3) 

k= 1 

In  Eq.  (3),  the  first  term  on  the  right-hand  side  estimates 
the  error  in  the  boundary  conditions  as  an  actual  value 
(BC)  minus  a  modeled  first  guess  (BC6).  The  second 
term  computes  the  weighted  residual  between  the  model 
and  the  observations,  a  sum  over  individual  innovations 
(yk  -  H^N)  at  each  location  k  (e.g.,  Kalnay  2003;  Park 
and  Xu  2009;  Lahoz  et  al.  2010;  Swinbank  et  al.  2003). 
In  this  term,  Hk  is  an  observational  operator,  which 
projects  the  model  solution  onto  the  location  of 


observation  y&.  Terms  QB^  and  R-1  are  weights  based 
on  the  relative  variance  of  each  error  source.  The  cost 
function  is  minimized  by  finding  the  set  of  controls 
(boundary  conditions)  for  which  the  first  variation  (or 
gradient)  of  J(N)  is  zero.  An  adjoint  variable  is  in¬ 
troduced  to  facilitate  the  computation  of  the  gradient  of 
the  cost  function  with  respect  to  the  controls.  This  leads 
to  a  set  of  coupled  Euler-Lagrange  equations  that  are 
solved  iteratively  (see  appendix  A  for  additional  de¬ 
tails).  For  the  present  analysis,  only  the  boundary  con¬ 
ditions  are  controlled,  which  for  the  steady-state  system  is 
equivalent  to  a  strong-constraint  approach  to  data  as¬ 
similation,  wherein  the  full  model  is  assumed  to  generate 
no  errors — that  is,  all  errors  are  propagated  into  the  do¬ 
main  from  the  boundary. 

The  assimilation  system  utilized  here  consists  of  the 
adjoint  to  linear  homogeneous  SWAN  coupled  to  an 
internal  copy  of  the  full  forward  SWAN  and  supple¬ 
mented  by  several  additional  modules  that  compute  and 
track  the  cost  function  and  optimize  system  convergence 
(Fig.  1).  The  system  is  wholly  separate  from  the  original 
forward  SWAN,  interfacing  only  at  the  beginning  and 
end  of  the  data  assimilation  process.  To  begin,  the  ex¬ 
ternal  forward  SWAN  is  initialized  with  boundary  data 
from  regional  models  and  used  to  compute  estimated 
wave  spectra  at  all  measurement  locations  in  the  do¬ 
main.  An  initial  estimate  for  the  cost  function  J(N)  is 
then  obtained  from  the  preliminary  boundary  condi¬ 
tions,  SWAN  spectra,  and  observed  spectra,  and  its 
gradient  is  calculated  at  the  offshore  boundary  using  the 
Euler-Lagrange  equations.  A  conjugate  gradient  tech¬ 
nique  (Polak  and  Ribiere  1969;  Walker  2006)  is  used  to 
adjust  the  estimate  toward  its  optimal  value.  If  the  cost 
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function  is  not  yet  minimized  by  the  adjusted  estimate, 
then  the  adjoint  model  is  solved  to  obtain  additional 
corrections  to  the  boundary  conditions,  and  the  internal 
forward  SWAN  is  run  with  the  new  boundary  inputs  to 
generate  revised  spectral  estimates.  The  revised  spectra 
are  used  to  recompute  the  cost  function,  and  the  itera¬ 
tions  continue  in  this  manner  until  J(N )  is  minimized 
and  the  optimal  corrected  boundary  conditions  are  ob¬ 
tained.  These  are  provided  to  the  external  forward 
SWAN,  which  computes  the  corrected  solution  for  the 
entire  domain. 

Adjoint  construction 

Individual  adjoint  SWAN  subroutines  are  built  with 
the  Parametric  FORTRAN  compiler  (PFC)  utility 
(Erwig  et  al.  2007),  which  uses  automatic  differentia¬ 
tion,  applying  standard  rules  in  combination  with  user- 
specified  active  and  dependent  variables  to  linearize 
or  transpose  each  forward  subroutine.  Active  and  de¬ 
pendent  variables  are  identified  in  the  forward  sub¬ 
routine  and  listed  in  a  parameter  file  that  is  read  by 
PFC.  For  example,  if  vector  c  is  active  in  the  iterative 
computation, 

do  i  =  2  :  nlocs  —  1 
a(7  +  1)  =  c (/  -  1)  +  B  o  c (/)  +  E 
enddo ,  (4) 

where  nlocs  is  the  number  of  locations  on  a  geographic 
grid,  a  is  a  dependent  variable,  B  is  a  matrix  of  co¬ 
efficients,  and  E  is  an  independent  fixed  constant,  then 
the  (perturbation  type)  tangent  linear  version  of  Eq.  (4) 
has  essentially  the  same  form,  except  that  the  unperturbed 
constant  E  is  missing  as  shown: 

do  i  =  2  :  nlocs  -  1 
Sa(/  +  1)  =  8c(i  —  1)  +  B°Sc(/) 
enddo .  (5) 

The  adjoint  of  Eq.  (4)  [and  (5)]  provides  the  sensitivity 
of  the  observation  location  with  respect  to  the  boundary 
control.  It  propagates  a  first  derivative  at  a  given  ob¬ 
servation  location  back  to  the  boundary  control,  taking 
into  account  the  model  state  and  physics.  It  can  be  used 
to  determine  how  the  error  in  an  estimate  of  c  at  a  given 
grid  location  is  related  to  a  corresponding  error  in  c  at 
the  boundary,  in  the  presence  of  the  field  specified  by 
Eq.  (4).  In  the  adjoint,  time  and  spatial  loops  are  re¬ 
versed  along  with  the  relationships  of  the  active  vari¬ 
ables,  the  matrix  of  coefficients  is  transposed,  and  the 
resulting  expression  for  each  adjoint  variable  becomes 
cumulative  over  the  loop  as  shown: 


do  i  =  (nlocs  -  1)  :  2 
ad_c(/  -  1)  =  ad_c(/  -  1)  +  ad_a(/  +  1) 
ad_c(/)  =  ad_c(/)  +  BT°ad_a(z'  +  1) 
enddo .  (6) 

In  the  adjoint  model  as  a  whole,  subroutines  themselves 
are  also  called  in  reverse  order  of  those  in  the  forward 
model.  Variables  that  are  not  active  or  dependent  are 
computed  in  the  adjoint  as  they  are  in  the  forward 
model,  and  they  must  have  the  same  values  as  their 
forward  counterparts  at  the  same  times  and  locations. 
Equation  (2)  and  its  adjoint  are  presented  in  discrete 
form  in  appendix  B. 

Subroutines  utilized  by  stationary  homogeneous 
SWAN  are  generally  already  linear,  so  the  tangent  lin¬ 
ear  (TL)  version  of  these  subroutines  is  essentially  un¬ 
changed  from  the  original  forward  version  [similar  to 
Eqs.  (4)  and  (5)].  Active  and  dependent  variables  in 
SWAN  normally  include  the  wave  action  and  anything 
that  depends  on  it.  Examples  of  independent  variables 
include  the  energy  propagation  velocity  C  and  the  water 
depth.  As  each  numerical  adjoint  subroutine  in  this 
system  is  derived  directly  from  its  forward  counterpart, 
the  principal  physics  of  the  original  subroutine  are  pre¬ 
served  and  the  adjoint  model  retains  the  associated 
properties  of  the  discretized  linear  forward  model.  To 
improve  overall  accuracy,  all  real  arrays  in  the  assimi¬ 
lation  system’s  adjoint  and  forward  SWAN  programs 
are  compiled  with  double  precision.  Using  these  methods, 
a  full  set  of  adjoint  subroutines  is  constructed  for  sta¬ 
tionary  homogeneous  SWAN. 

3.  Validation 

In  matrix  notation,  the  action  of  a  given  linearized 
SWAN  subroutine  upon  an  input  array  u ,  generating  an 
output  array  v,  may  be  expressed  as  Au  =  v,  while  that  of 
the  corresponding  adjoint  subroutine  can  be  written  as 
ATv  =  u.  For  a  real  matrix  A  with  transpose  AT  and  ar¬ 
bitrary  real  vectors  u  and  v,  the  following  inner  product 
is  an  identity: 

(Au,v)  =  (u,  Atv).  (7) 

For  a  given  subroutine  and  its  adjoint,  this  implies  that  if 
we  initialize  the  forward  subroutine  with  an  arbitrary 
vector  u  (obtaining  an  output  Au)  and  then  initialize  its 
adjoint  with  a  second  arbitrary  vector  v  (obtaining  an 
output  ATv),  then  an  inner  product  of  the  forward  rou¬ 
tine  output  with  the  adjoint  input  should  have  the  same 
value  as  an  inner  product  of  the  adjoint  output  with  the 
forward  input. 
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This  identity  is  also  exhibited  in  the  symmetry  of  the 
representer  matrix  (Bennett  2002),  which  provides 
a  convenient  method  for  testing  the  consistency  of  the 
complete  adjoint  program  and  its  forward  counterpart. 
For  a  given  number  m  of  observation  locations  (in  the 
spatial-spectral  domain),  the  representer  function  for 
each  location  is  obtained  by  initializing  the  adjoint  with 
a  unit  impulse  centered  at  that  location.  The  adjoint 
solution  (at  all  grid  points)  is  then  used  to  force  the 
forward  model  (running  in  stationary  mode  with  no 
sources  or  sinks).  The  forward  SWAN  solutions  at  the  m 
observation  locations  then  fill  one  column  of  the  m-by-m 
representer  matrix  (HAAtHt,  where  H  is  composed  of 
m  observational  operators,  Hk ,  for  k  =  1  to  m).  If  the 
adjoint  is  consistent  with  the  stationary  linear  forward 
model,  then  the  representer  matrix  should  be  symmetric 
to  machine  precision  (Ngodock  and  Carrier  2013). 

Subroutines  in  the  new  adjoint  are  validated  individ¬ 
ually  by  means  of  the  above-mentioned  inner  product 
test  [Eq.  (7)].  In  each  test,  input  vectors  u  and  v  are  as¬ 
signed  different  sets  of  random  values,  and  inner  prod¬ 
uct  values  are  checked  at  multiple  locations  on  an 
idealized  bathymetry.  For  each  individual  subroutine  or 
group  of  related  subroutines  tested,  Eq.  (7)  is  consis¬ 
tently  satisfied  to  within  machine  precision.  Specifically, 
for  all  such  tests, 

(Am..)  — (u.Atv)s5x10-i, 

(Au,  v) 

Consistency  of  the  full  adjoint  code  is  checked  with  an 
impulse  test,  in  which  the  representer  matrix  is  com¬ 
puted  for  a  limited  number  of  points,  as  outlined  above, 
and  its  symmetry  is  checked.  The  tests  are  conducted 
using  bathymetry  from  Santa  Rosa  Island,  Florida  (Fig.  2), 
and  selecting  four  grid  locations  (TA1,  TA2,  SAB,  SIB) 
matching  those  of  actual  instruments  in  a  2009  experiment 
(Edwards  et  al.  2009;  Veeramony  et  al.  2010).  Locations 
TA1  and  SAB  are  closer  to  the  shoreline,  where  bathy¬ 
metry  is  alongshore  uniform  and  some  wave  shoaling 
and  breaking  may  occur,  while  TA2  and  SIB  are  outside 
the  surf  zone,  near  the  southwestern  and  southeastern 
grid  boundaries,  respectively.  In  four  separate  subtests, 
the  adjoint  is  initialized  with  a  unit  impulse  at  a  single 
frequency  and  direction  for  each  of  the  four  locations 
described  above.  For  example,  at  location  TA1  the 
initialization  is 

v(x tar  yTAVf  =  0.066Hz,  0  =  50°)  =  1,  v=0  elsewhere. 

In  each  subtest,  the  adjoint  is  allowed  to  run  for  a  single 
iteration.  Its  output  is  used  to  initialize  forward  SWAN 


Fig.  2.  Model  bathymetry  for  inner  product,  impulse  and  twin 
experiment  tests  based  on  measurements  at  Santa  Rosa  Island,  FL, 
recorded  in  2009.  Direction  on  the  y  axis  is  north.  Shoreline  is 
shaded.  The  four  grid  locations  used  for  impulse  and  twin  tests  are 
circled  and  labeled. 

(in  stationary  linear  mode),  whose  output  is  in  turn 
recorded  at  all  four  of  the  selected  locations.  The  results 
are  compiled  into  a  (4  X  4)  subsection  of  the  representer 
matrix,  which  is  found  to  be  symmetric  at  all  locations 
within  model  accuracy  of  order  lO-12,  thus  confirming 
that  the  full  adjoint  is  consistent  with  the  stationary 
linear  forward  model. 

4.  Simulations 

Two  sets  of  simulations  are  conducted  to  further  ex¬ 
amine  the  performance  of  the  numerical  adjoint.  First, 
the  numerical  adjoint  approach  as  a  whole  is  evaluated 
by  direct  comparison  with  the  Walker  (2006)  approach, 
using  a  twin  experiment  where  artificial  observations 
are  generated  by  the  full  forward  model  with  idealized 
boundary  conditions.  Following  this,  the  assimilation 
system  is  evaluated  with  measured  spectra  from  multiple 
instruments  at  Duck,  North  Carolina.  For  both  tests, 
errors  in  boundary  control  are  assumed  to  be  completely 
spatially  correlated  and  uniform.  Final  spectral  esti¬ 
mates  are  compared  to  observations  at  boundary  and 
interior  grid  locations. 

a.  Twin  experiment 

The  idealized  twin  experiment  validation  tests  are 
based  directly  on  those  conducted  by  Veeramony  et  al. 
(2010)  as  part  of  their  original  evaluation  of  the  Walker 
(2006)  adjoint,  utilizing  data  and  bathymetry  from  the 
aforementioned  2009  experiment  at  Santa  Rosa  Island, 
Florida.  The  fully  nonlinear  forward  SWAN,  initialized 
at  seaward  and  lateral  boundaries  with  a  parametric  Joint 
North  Sea  Wave  Project  (JONSWAP)  wave  spectrum,  is 
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Fig.  3.  Evolution  of  (left)  cost  function  and  (right)  norm  of  gradient  of  cost  function  plotted  vs  assimilation  iteration  for  all  twin 
experiment  tests.  Asterisks  and  circles  represent  numeric  and  analytic  adjoint  system  results,  respectively.  All  values  are  normalized  by 
the  value  at  the  first  iteration.  Heavy  dashed  line  (left)  marks  the  0.75%  criterion  for  the  cost  function  below  which  each  assimilation  is 
stopped.  Note  log  format  on  the  y  axis. 


used  to  compute  pseudo-observed  spectra  at  the  same  four 
locations  used  in  the  2010  analysis  (Fig.  2).  JONSWAP 
parameters  are  based  on  prevailing  waves  during  the 
experiment,  which  were  relatively  mild — significant 
wave  height  (Hs)  <1.5  m — and  approached  generally 
from  the  southeast  and  normal  to  the  offshore  grid 
boundary.  Four  of  the  six  cases  from  Veeramony  et  al. 
(2010)  are  rerun  for  both  assimilation  systems  (cases  4 
and  6,  featuring  local  wind  forcing  nonlinearities,  are 
excluded).  For  each  case,  four  separate  assimilation  tests 
are  conducted  with  each  adjoint  system.  In  the  first  of 
these  tests,  the  adjoint  is  initialized  with  the  innovation 
at  grid  location  TA1;  in  the  second  test,  at  TA2;  in  the 
third  test,  at  SAB;  and  in  the  fourth  test,  at  SIB.  The 
system  assumes  a  zero-energy  state  as  a  first  guess  (i.e., 
spectra  equal  to  zero  at  boundary  and,  consequently, 
throughout  the  grid),  so  the  initializing  innovation  for  a 
given  test  (i.e.,  the  observed  minus  the  estimated  spec¬ 
trum)  is  simply  equal  to  the  observed  spectrum  at  that 
location. 

In  each  test,  adjoint-computed  values  of  spectral 
density  at  the  boundary  nodes  are  averaged  across  all 
locations  on  the  seaward  boundary  and  provided  as 
uniform  offshore  boundary  conditions  to  the  assimila¬ 
tion  system’s  internal  forward  SWAN  [with  wave 
breaking  and  other  source/sink  subroutines  activated, 
duplicating  the  setup  of  Veeramony  et  al.  (2010)  for 
each  case].  Along  lateral  boundaries,  spectra  are  set  to 
a  tapering  fraction  of  the  full  offshore  spectrum,  ranging 
from  100%  at  the  seaward  end  to  zero  at  the  shoreward 
end  (SWAN  Team  2011).  A  revised  estimate  for  the 


wave  spectrum  at  the  initializing  location  is  obtained 
and  compared  with  the  original  data  spectrum.  If  the 
cost  function  is  not  yet  minimized,  then  the  adjoint  is 
reinitialized  with  a  recomputed  innovation.  Each  as¬ 
similation  system  is  allowed  to  run  until  its  cost  function 
declines  to  0.75%  of  its  original  value.  By  this  point  the 
associated  cost  function  gradients  have  been  reduced  to 
less  than  5%  of  their  initial  values  and  are  generally 
flattening  out,  as  are  the  cost  function  values  themselves 
(Fig.  3). 

In  each  test,  the  two  assimilation  systems’  relative 
performance  is  qualitatively  evaluated  by  comparing  how 
well  they  reproduce  the  observed  spectra  at  the  offshore 
boundary  and  all  four  instrument  locations,  including 
nonassimilated  spectra  as  well  as  the  selected  innovation 
spectrum.  For  a  more  quantitative  comparison,  overall 
model  accuracy  is  also  evaluated  using  an  RMS  skill  score 
computed  from  spectral  densities  as  shown: 


l[Smod(U°j)-Sobs(U°j)f 

- -  - .  (8) 

Jl[sobs(U°j)f 

V  ij 

In  Eq.  (8),  Smod  is  the  model  spectrum  and  Sobs  is  the 
observed  spectrum  (from  nonlinear  forward  SWAN). 
Spectral  energy  densities  are  first  squared  and  then 
summed  over  all  frequencies  (ft)  and  directions  (0y). 

In  general,  spectral  estimates  from  the  present  assimi¬ 
lation  system  are  very  similar  to  those  from  the  Walker 
(2006)  system.  Both  systems  have  errors  of  similar 


skill  =  1  - 
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Fig.  4.  (top  panels)  Twin  experiment  results,  comparing  stage  1  numerical  adjoint  to  Walker  (2006)  analytical  adjoint.  Both  are  ini¬ 
tialized  with  data  from  nonlinear  forward  SWAN  at  the  TAl  location.  Four-panel  plots  for  each  of  the  four  interior  locations  (TAl,  TA2, 
SAB,  and  SIB)  include  spectral  results  from  numeric/analytic-adjoint-based  assimilation  systems  (top  two  panels  of  each  set)  and  the 
difference  between  those  results  and  the  original  SWAN  spectra  (bottom  two  panels).  The  three  panels  below  the  grid  map  show  mean 
predicted  boundary  spectrum:  numeric,  analytic,  and  numeric  -  analytic,  (bottom)  Comparison  of  normalized  distribution  of  adjoint 
boundary  energy  along  offshore  boundary  for  the  two  adjoint  types. 


magnitude  when  compared  to  observed  spectra  (e.g., 
Fig.  4).  When  skill  values  are  computed  and  averaged 
over  all  twin  experiment  tests  (64  computations  per 
adjoint),  the  numerical  assimilation  system  achieves 
a  mean  skill  of  0.914  (range  =  0.82-0.98),  roughly 
equivalent  to  the  analytical  adjoint  system’s  mean  of 
0.911  (range  =  0.78-0.97). 

b.  Observed  data — Duck,  North  Carolina 

The  final  evaluation  of  the  simplified  numerical  ad¬ 
joint  tests  how  well  the  assimilation  system  can  re¬ 
capture  actual  measured  spectra,  starting  from  an 
essentially  arbitrary  first  guess.  Spectral  wave  datasets 
are  obtained  from  the  Army  Corps  of  Engineers’  Field 
Research  Facility  (FRF)  in  Duck,  North  Carolina,  for 
several  periods  that  also  include  nearshore  bathymetry 
measurements.  To  minimize  nonlinear  effects,  study 
dates  are  selected  from  relatively  mild  wave  climates  in 


which  significant  wave  heights  are  less  than  1.5  m  (Fig.  5). 
Bathymetry  data  recorded  for  the  Duck  “minigrid”  on 
16  April,  1  June,  and  29  July  2010  are  linearly  extrapo¬ 
lated  alongshore  to  5  times  the  minigrid  width  and  off¬ 
shore  to  17-m  depth.  Spectral  data  are  available  from  up 
to  six  nearshore  instruments  (Fig.  6),  including  four 
Nortek  acoustic  wave  and  current  (AWAC)  profilers  at 
depths  of  5,  6,  8,  and  11  m,  and  the  FRF’s  8-m  array  of 
15  near-bottom-mounted  pressure  sensors.  Boundary 
spectra  for  the  model  grid  are  provided  by  a  Datawell 
Waverider  630  buoy,  moored  approximately  3  km  off¬ 
shore  in  roughly  17-m  water  depth.  For  the  16  April  and 
1  June  simulations,  four  instruments  are  used  (the  5  and 
8  m  AWAC  profilers  were  not  available).  For  29  July,  all 
six  instruments  are  used. 

Owing  to  the  manner  in  which  they  are  derived, 
spectra  from  the  four  AWACs  have  lower  resolution 
than  those  from  the  8-m  array  and  the  Waverider,  and 
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16  Apr  1  Jun  29Jul 


Fig.  5.  (left  to  right)  Time  series  of  (top  to  bottom)  significant  wave  height,  peak  period,  and  peak  direction  as  measured  at  the  3-km 
Waverider  buoy  for  periods  roughly  concurrent  with  the  three  bathymetry  measurements  included  in  these  tests.  Conditions  at  times  of 
selected  initializing  spectra  are  marked  with  asterisks  in  each  panel. 


they  also  tend  to  have  a  significantly  broader  directional 
spread,  with  a  small-to-moderate  amount  of  wave  en¬ 
ergy  directed  offshore  at  some  sensors.  It  is  not  known 
how  much  of  this  additional  spreading  is  a  result  of  in¬ 
strument  error  and  how  much  is  due  to  unrepresented 
nonlinear  effects  such  as  wave  breaking  and  reflection. 
The  current  assimilation  system  does  not  presently  as¬ 
similate  offshore-directed  wave  energy.  To  more  effec¬ 
tively  evaluate  actual  system  capabilities,  assimilated 
spectra  are  therefore  limited  to  onshore  directions  in 
these  tests  (i.e.,  all  offshore-directed  spectral  energy — 
generally  much  less  than  one  percent  of  total  energy — is 
set  to  zero  in  observed  spectra).  Even  with  this  limita¬ 
tion,  the  varied  directional  spread  among  initializing 
instruments  has  a  noticeable  effect  on  model  output  that 
will  be  discussed  further  below. 

The  present  tests  are  designed  to  evaluate  the  system’s 
ability  to  assimilate  real  spectral  data  and  investigate  its 


sensitivity  to  the  observed  data  locations.  For  each  date, 
the  adjoint  is  initialized  with  innovations  at  one  or  more 
locations  in  nine  different  cases.  In  case  1,  only  spectra 
from  the  FRF  8-m  array  are  used.  In  case  2,  the  adjoint 
is  initialized  at  all  available  AWAC  profiler  locations, 
with  each  innovation  equally  weighted.  Case  3  utilizes 
equally  weighted  spectra  from  all  (three  or  five)  interior 
grid  locations.  In  case  4,  all  interior  spectra  are  again 
used,  but  the  8-m  array  is  given  additional  weight.  Case  5 
uses  all  interior  spectra  and  instead  gives  additional 
weight  to  each  of  the  AWAC  profilers  (distributed 
evenly  among  them).  Cases  6,  7,  8,  and  9  are  initialized 
only  at  the  5,  6,  8,  and  11  m  AWAC  profilers,  respec¬ 
tively  (where  available). 

As  with  the  twin  experiment,  a  zero-energy  spectrum 
is  again  used  for  the  first  guess,  so  that  at  each  location 
the  innovation  is  just  equal  to  the  observed  spectrum  for 
the  first  iteration.  Wave  breaking  is  activated  for  the 
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Fig.  6.  Map  of  nearshore  instrument  locations  at  FRF  (offshore  Waverider  buoys  not  shown).  Instruments  used  in  this  study  are  circled. 

(Image  courtesy  of  FRF,  http://www.frf.usace.army.mil/frfzoom.shtml.) 


system’s  internal  forward  SWAN  model  to  maintain 
system  stability,  but  other  nonlinear  sources  and  sinks 
are  not  activated.  Each  simulation  is  allowed  to  run  until 
formal  convergence  has  been  achieved  (i.e.,  the  cost 
function  value  has  become  essentially  constant).  Final 
estimated  spectra  are  recorded  at  the  seaward  boundary 
and  at  all  observation  locations,  and  wave  statistics 


including  significant  wave  height,  mean  direction,  mean 
period,  and  directional  spread  are  computed  from  each 
modeled  and  observed  spectrum. 

Model  skill  [evaluated  using  Eq.  (8)]  ranges  from 
a  low  near  zero  to  a  high  of  0.92,  with  a  mean  skill  of  0.52. 
Model  results  for  the  test  series  are  summarized  in  Table  1. 
The  highest  skill  values  (averaging  0.73)  are  obtained  for 


Table  1.  Adjoint  skill  scores  for  tests  at  Duck,  NC.  Tests  were  conducted  for  three  separate  dates  in  2010  (16  Apr,  1  Jun,  and  29  Jul  7). 
Where  data  were  available  from  multiple  dates,  a  range  of  skill  scores  is  provided  in  the  table.  First  column  shows  case  number  and  second 
column  shows  initializing  location(s)  for  the  adjoint  system.  AWACs  implies  all  available  AWAC  sensors  were  used;  All  implies  all 
available  sensors  were  used.  Overwt  is  used  for  cases  where  one  or  more  sensors  were  overweighted  relative  to  other  sensors  used. 


Case 

Initializing  location 

Output  location 

FRF 

AWAC  (5  m) 

AWAC  (6  m) 

AWAC  (8  m) 

AWAC  (11  m) 

1 

FRF 

0.67-0.92 

0.42 

0.19-0.50 

0.26 

0.00-0.51 

2 

AWACs 

0.32-0.61 

0.53 

0.61-0.67 

0.52 

0.63-0.67 

3 

All 

0.50-0.73 

0.55 

0.60-0.67 

0.48 

0.48-0.64 

4 

All  Overwt  FRF 

0.76-0.91 

0.44 

0.32-0.53 

0.15 

0.26 

5 

All  Overwt  AWACs 

0.37-0.65 

0.55 

0.63-0.67 

0.51 

0.08-0.62 

6 

AWAC  (5  m) 

0.47 

0.71 

0.52 

0.29 

0.27 

7 

AWAC  (6  m) 

0.35-0.58 

0.51 

0.68-0.72 

0.44 

0.45-0.56 

8 

AWAC  (8  m) 

0.39 

0.45 

0.56 

0.58 

0.60 

9 

AWAC  (11  m) 

0.26-0.53 

0.39 

0.44-0.55 

0.47 

0.70-0.79 
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Fig.  7.  Evolution  of  (left)  cost  function  and  (right)  norm  of  gradient  of  cost  function,  plotted  vs  assimilation  iteration  for  all  adjoint 
simulations,  with  colors  and  symbols  representing  the  nine  different  cases  shown  in  Table  1.  All  values  are  normalized  by  the  value  at  the 
first  iteration.  Note  log  format  on  y  axis. 


tests  in  which  a  single  location  (e.g.,  the  8-m  array)  is 
used  for  both  input  and  output.  In  general,  spectral  es¬ 
timates  are  poor  for  cases  where  the  initializing  in¬ 
strument  type  differs  from  the  instrument  type  at  the 
output  location.  The  worst  results  (average  skill  =  0.38) 
are  obtained  for  spectral  estimates  at  AWAC  profiler 
locations  in  cases  1  and  4,  when  the  model  is  initialized 
either  solely  at  the  8-m  array  or  at  all  locations  with  the 
8-m  array  overweighted.  The  relative  change  in  the  cost 
function  and  its  gradient  is  shown  for  all  cases  in  Fig.  7. 
Cost  functions  and  gradients  decrease  most  rapidly  for 
tests  initialized  at  the  farthest  offshore  sensors  (cases  1 
and  9),  particularly  when  waves  are  small  (16  April). 
The  slowest  convergence  is  obtained  when  using  data 
from  just  the  shallowest  AWAC  profiler  (case  6)  on  the 
largest  wave  day  (29  July),  while  the  smallest  reduction 
of  the  cost  function  is  obtained  for  equally  weighted 
tests  initialized  at  all  available  sensors  (case  3). 

Modeled  and  observed  spectral  statistics  are  compared 
at  all  grid  locations  for  all  cases  in  Fig.  8.  Regardless  of 
initializing  data  and  weights,  the  model  consistently 
underpredicts  significant  wave  height.  The  assimilation 
system  appears  to  perform  slightly  better  for  cases  with 
multiple  initializing  sensors.  This  is  not  surprising;  using 
a  larger  number  of  initializing  data  locations  of  compa¬ 
rable  quality  will  generally  improve  modeled  wave 
heights.  Although  they  have  lower  spectral  resolution, 
the  AWAC  profilers  capture  spectral  moments  roughly 
as  well  as  the  8-m  array  (K.  Hathaway  2011,  personal 
communication).  Wave  height  results  appear  to  be  some¬ 
what  better  for  smaller  waves,  suggesting  that  neglected 
effects  such  as  wind  forcing,  wave  breaking,  and  triad  in¬ 
teractions  (per  their  formulations  in  the  SWAN  model) 
likely  play  a  role  in  generating  these  errors. 


Mean  wave  directions  are  recovered  fairly  well,  al¬ 
though  model  estimates  are  slightly  closer  to  shore 
normal  (90°)  than  observations.  Mean  periods  are  cap¬ 
tured  quite  well  in  nearly  all  cases,  with  a  very  high 
correlation  and  consistent  1:1  ratio  of  estimates  to  ob¬ 
servations.  Some  error  may  result  from  neglected  wind 
forcing  and  triad/quadruplet  interactions,  which  shift 
observed  wave  energy  to  different  frequencies  (with 
different  refractive  properties).  Directional  spreading 
results  are  poor,  largely  owing  to  differences  in  the 
spectra  produced  by  each  type  of  instrument.  Predicted 
directional  spreads  at  AWAC  profiler  locations  in  cases 
1  and  4  are  consistently  lower  than  observed  values  at 
the  AWAC  profilers.  This  is  expected  because  of  the 
stronger  effects  of  the  directionally  narrower  8-m  array 
spectra  on  model  output  in  these  cases.  Inversely,  in 
cases  2  and  5-9,  modeled  directional  spreads  are  over¬ 
estimated  at  the  8-m  array  location  because  the  adjoint 
is  more  strongly  influenced  by  the  directionally  broad 
AWAC  profilers. 

5.  Discussion 

Twin  experiment  spectral  estimates  from  the  data  as¬ 
similation  system  that  uses  the  numerical  adjoint  are 
consistently  as  good  as  those  from  the  Walker  (2006) 
system  built  around  an  analytical  adjoint,  and  the  two 
models  show  equal  levels  of  skill.  Results  from  the  nu¬ 
merical  adjoint  tests  with  observed  data  at  Duck,  North 
Carolina,  are  good  for  small  wave  cases  and  become  fair 
or  worse  as  spectral  energies  increase.  These  comparisons 
illustrate  some  of  the  shortcomings  of  using  the  reduced 
forward  model  to  derive  the  adjoint,  particularly  in  esti¬ 
mates  of  significant  wave  height  and  directional  spread. 
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Fig.  8.  Summary  plots  of  modeled  vs  observed  spectral  statistics  for  all  cases,  showing  results  at  interior  and  boundary  data  locations: 
(top)  (left)  Hs  and  (right)  mean  direction;  and  (bottom  (left)  mean  period  and  (right)  direction  spread.  Nine  symbols  represent  each  of  the 
initializing  scenarios  described  in  the  text.  The  solid  line  in  each  panel  represents  perfect  agreement  between  model  and  data.  The  dashed 
line  shows  the  best  linear  fit  to  the  data;  its  r2  value  is  provided  in  each  panel.  For  the  mean  direction  panel,  shore  normal  is  90°. 


Nonlinear  sources  of  wave  energy  from  wind  forcing 
and  wave-wave  interactions  are  not  represented  in  the 
adjoint  code,  although  using  the  numerical  approach 
to  derive  the  adjoint  will  ultimately  make  this  repre¬ 
sentation  possible.  Wave  breaking  was  activated  in  the 
system’s  internal  forward  SWAN  model,  but  its  coun¬ 
terpart  was  not  available  in  the  adjoint.  This  mismatch 
between  the  forward  and  adjoint  SWAN  components  of 
this  suboptimal  assimilation  system  may  contribute  to  the 
underestimation  of  larger  wave  heights  by  the  model. 
Directional  spreading  also  appears  to  increase  at  the 
shallower  AWAC  profiler  sensors,  at  least  partially  owing 
to  triad  and  quadruplet  interactions  in  or  near  the  surf 
zone,  both  of  which  are  also  neglected  in  the  numerical 
adjoint. 

To  a  considerable  extent,  however,  system  perfor¬ 
mance  is  hampered  by  the  use  of  spectra  from  different 


types  of  instruments.  Spectral  processing  differences  be¬ 
tween  the  AWAC  profilers  and  the  8-m  array  significantly 
limit  the  adjoint’s  performance  in  matching  observed 
spectra  (e.g.,  Fig.  9,  left  and  center  panels).  Adjoint  esti¬ 
mates  in  these  tests  might  be  improved  if  the  directional 
spreads  of  all  input  spectra  could  be  normalized  (e.g., 
using  the  relative  directional  spreads  measured  concur¬ 
rently  at  the  8-m  array  and  the  8-m  AWAC  profiler). 
However,  this  modification  would  have  to  be  based  on  an 
explicit  mathematical  representation  of  the  complex  re¬ 
lationship  between  the  different  techniques  used  to  derive 
spectra  at  each  instrument,  which  may  or  may  not  be 
derivable.  Spectra  from  the  FRF  8-m  array  of  15  pressure 
sensors  are  obtained  with  an  iterative  maximum  likeli¬ 
hood  estimator  (IMLE)  applied  to  8192-s  time  series 
(US ACE  Field  Research  Facility  2011).  In  contrast, 
AWAC  profiler  spectra  are  generated  using  the  combined 
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Fig.  9.  (left)  Observed  and  (middle)  modeled  frequency-directional  spectra  at  the  (top)  FRF 
array  and  (bottom)  8-m  AWAC  profiler  for  29  Jul  2010.  (right)  Corresponding  modeled  and 
observed  frequency  spectra  for  each  instrument  (log  scale  ony  axis).  Skill  values  for  2D  and  ID 
estimates  are  provided  in  parentheses.  Model  results  are  for  case  1,  in  which  the  adjoint  was 
initialized  only  with  a  spectrum  from  the  FRF  8-m  array. 


three-beam  UV  velocity  data  together  with  pressure 
measurements  in  another  IMLE  analysis,  with  slightly 
different  methods  used  for  low-  and  high-frequency 
ranges.  The  offshore  Waverider  buoy  employs  Fourier 
transform  techniques  with  time  series  of  vertical  and 
horizontal  buoy  displacements,  and  an  MLE  analysis 
is  applied  to  the  resulting  Fourier  coefficients  to  ob¬ 
tain  directional  spectra  (K.  Hathaway  2011,  personal 
communication) . 

To  remove  directional  spreading  from  consideration 
in  evaluating  model  performance,  a  second  set  of  skill 
calculations  is  made  using  modeled  and  observed  fre¬ 
quency  spectra  instead  of  frequency-directional  spectra. 
The  revised  skill  score  expression  is 


lisjfd-sM)]2 

'  /v  (9) 

yi  ysQm2 

where  Smo d  and  Sohs  from  Eq.  (8)  have  been  summed 
over  directions  to  obtain  Sm  and  S0  above,  respectively. 
The  resulting  skill  scores  (Table  2)  are  significantly 


better  than  those  obtained  with  2D  spectra.  For  calcu¬ 
lations  based  on  Eq.  (9),  model  skill  at  Duck,  North 
Carolina,  now  ranges  from  0.51  to  0.95  with  a  mean 
overall  skill  of  0.71.  This  result  suggests  that,  although 
directional  spreading  differences  are  dominant,  there 
are  also  some  frequency  domain  variations  between 
spectra  from  the  8-m  array  and  the  AWAC  profilers.  For 
the  sample  case  in  Fig.  9,  a  comparison  of  ID  frequency 
spectra  (right  two  panels)  illustrates  the  considerable 
improvement  in  model  performance  when  directional 
spread  is  not  considered.  In  this  case,  while  model  skill  at 
the  8-m  array  improves  marginally  (from  0.88  to  0.89), 
there  is  a  much  more  substantial  increase  at  the  8-m 
AWAC  profiler  (from  0.26  to  0.79).  Remaining  spectral 
differences  at  the  8-m  AWAC  profiler  are  likely  partly 
due  to  instrumental  variations  and  in  part  to  nonlinear 
effects  that  are  neglected  in  the  adjoint. 

To  further  examine  the  importance  of  instrument 
type,  the  tests  in  section  4b  are  rerun,  replacing  observed 
data  at  each  location  with  artificial  spectra  generated 
using  fully  nonlinear  forward  SWAN.  The  forward 
model  is  initialized  with  data  from  the  offshore  Waver¬ 
ider  buoy,  and  wave  breaking,  triad  interactions,  wind 
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Table  2.  As  in  Table  1,  but  for  adjoint  skill  scores  for  tests  at  Duck,  NC,  computed  in  frequency  space  only  using  Eq.  (9). 


Output  location 


Case 

Initializing  location 

FRF 

AWAC  (5  m) 

AWAC  (6  m) 

AWAC  (8  m) 

AWAC  (11  m) 

1 

FRF 

0.83-0.92 

0.56 

0.59-0.67 

0.79 

0.51-0.77 

2 

AWACs 

0.68-0.76 

0.63 

0.67-0.84 

0.74 

0.67-0.74 

3 

All 

0.73-0.85 

0.66 

0.69-0.80 

0.78 

0.62-0.78 

4 

All  Overwt  FRF 

0.82-0.95 

0.70 

0.62-0.76 

0.85 

0.54-0.83 

5 

All  Overwt  AWACs 

0.71-0.78 

0.66 

0.70-0.83 

0.77 

0.66-0.77 

6 

AWAC  (5  m) 

0.67 

0.74 

0.69 

0.74 

0.72 

7 

AWAC  (6  m) 

0.57-0.83 

0.64 

0.74-0.76 

0.73 

0.55-0.71 

8 

AWAC  (8  m) 

0.63 

0.55 

0.60 

0.71 

0.66 

9 

AWAC  (11  m) 

0.51-0.78 

0.63 

0.66-0.76 

0.72 

0.78-0.83 

forcing  (10  m  s-1,  shore-normal),  and  quadruplet  in¬ 
teractions  are  all  activated.  Revised  assimilation  system 
estimates  of  wave  statistics  are  plotted  versus  those  from 
the  new  set  of  SWAN-based  “uniform  instrument”  ob¬ 
servations  in  Fig.  10.  The  results  are  significantly  better 
than  those  presented  in  Fig.  8,  especially  with  respect  to 
directional  spreading  estimates  (bottom-right  panel). 
Using  the  idealized  input  spectra,  model  skill  now 
ranges  from  0.54  to  0.93  with  a  mean  overall  skill  of  0.82 
(Table  3),  a  considerable  improvement  upon  the  levels 
seen  with  the  FRF  dataset  and  also  somewhat  better,  on 


Hs  (m) 


Mean  Period  (s) 


average,  than  the  results  obtained  with  ID  frequency 
spectra.  When  ID  skill  values  are  computed  for  these 
results  using  Eq.  (9),  model  skill  ranges  from  0.68  to 
0.93,  with  a  mean  skill  of  0.87. 

6.  Conclusions 

A  discrete  numerical  adjoint  to  the  wave  model 
SWAN  has  been  constructed  by  creating  individual 
numerical  adjoints  to  all  subroutines  in  the  stationary, 
homogeneous  part  of  the  forward  model.  The  adjoint 


Dir  Spread  (deg) 


Fig.  10.  As  in  Fig.  8,  but  for  summary  plots  of  modeled  vs  nonlinear  S WAN-generated  “observed”  spectral  statistics 
for  all  cases,  showing  results  at  interior  and  boundary  data  locations. 
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Table  3.  As  in  Table  1,  but  for  adjoint  skill  scores  for  tests  at  Duck,  NC,  using  “same  instrument”  pseudodata  observations  generated  by 

nonlinear  forward  SWAN  in  place  of  measured  data  at  each  location. 


Case 

Initializing  location 

Output  location 

FRF 

AWAC  (5  m) 

AWAC  (6  m) 

AWAC  (8  m) 

AWAC  (1  m) 

1 

FRF 

0.82-0.91 

0.93 

0.88-0.91 

0.89 

0.76-0.86 

2 

AWACs 

0.84-0.89 

0.92 

0.85-0.90 

0.88 

0.80-0.87 

3 

All 

0.84-0.90 

0.92 

0.85-0.90 

0.89 

0.79-0.85 

4 

All  Overwt  FRF 

0.84-0.89 

0.93 

0.84-0.91 

0.89 

0.77-0.86 

5 

All  Overwt  AWACs 

0.84-0.89 

0.92 

0.85-0.90 

0.88 

0.80-0.87 

6 

AWAC  (5  m) 

0.84 

0.90 

0.87 

0.84 

0.80 

7 

AWAC  (6  m) 

0.78-0.88 

0.92 

0.82-0.90 

0.87 

0.71-0.83 

8 

AWAC  (8  m) 

0.89 

0.93 

0.91 

0.89 

0.86 

9 

AWAC  (11) 

0.83-0.91 

0.93 

0.79-0.92 

0.91 

0.83-0.91 

and  its  subroutines  are  validated  with  standard  dot- 
product  and  impulse  tests,  after  which  the  adjoint  is  in¬ 
corporated  into  an  assimilation  system  together  with 
a  copy  of  the  original  forward  SWAN  code  and  modules 
for  minimizing  the  cost  function.  The  assimilation  sys¬ 
tem  is  then  tested  in  twin  experiments  with  artificial 
observed  spectra  (generated  by  nonlinear  SWAN)  at 
Santa  Rosa  Island,  Florida,  and  assimilations  with  actual 
measured  spectra  at  Duck,  North  Carolina.  In  the  twin 
experiments,  the  numerical  adjoint  performs  compa¬ 
rably  to  the  analytical  adjoint.  The  mean  skill  level 
[Eq.  (8)]  for  tests  with  artificial  data  is  about  0.91,  while 
the  average  skill  in  tests  with  measured  data  is  0.52. 
Differences  in  spectral  resolution  and  quality  at  the 
different  types  of  nearshore  instruments  in  Duck  play 
an  important  role  in  reducing  skill  levels  for  the  latter 
tests.  When  skill  scores  at  Duck  are  recomputed  for  ID 
frequency  spectra  only,  the  average  skill  value  in¬ 
creases  to  0.71.  Replacing  measured  spectra  in  the 
Duck  tests  with  artificial  “observations”  from  non¬ 
linear  SWAN  improves  the  system’s  average  2D  skill 
score  to  0.82  and  its  average  ID  skill  score  to  0.87.  The 
suboptimal  numerical  adjoint  system’s  performance  is 
sufficient  to  justify  the  development  of  a  complete  ad¬ 
joint  to  fully  nonlinear,  nonstationary  SWAN,  and  such 
an  effort  is  now  underway. 


APPENDIX  A 

Variational  Data  Assimilation 

The  following  material  is  included  for  readers  who  are 
unfamiliar  with  representer-based  data  assimilation  and 
follows  the  approach  of  Bennett  (2002).  Including  errors 
in  the  forcing,  the  initial  conditions,  and  the  dynamics, 
the  nonstationary,  inhomogeneous  SWAN  model  equa¬ 
tion  [Eq.  (1)]  may  be  written  in  one  representative  spatial 
dimension  as 


dN  dN 
—  +  c — 
dt  dx 


F+f 


u(x,  0)  =  I(x)  +  i(x) 


u(0,t)  =  B(t)  +  b(t), 


(Al) 


where  N  =  N(x,  t)  and  /;  i  and  b  represent  the  errors  in 
the  model  dynamics,  initial  and  boundary  conditions 
respectively,  with  covariances  Cf(x,  t,x'  ,t'),  Q(x,  x'),  and 
Cb(t,  t').  Here,  the  idealized  single  spatial  dimension  x 
represents  a  combination  of  two  spatial  dimensions 
(X,  Y)  and  two  additional  spectral  dimensions  (/,  6) 
from  the  “real  world.”  It  is  assumed  that,  during  the 
time  interval  [0,  7],  observations  of  part  of  the  state  are 
collected,  having  the  form 


Acknowledgments.  This  work  was  funded  by  the  Na¬ 
tional  Academy  of  Sciences  and  the  Naval  Research 
Laboratory,  through  a  National  Research  Council  Post¬ 
doctoral  Research  Associateship,  and  by  the  Office  of 
Naval  Research  through  the  6.2  NRL  Core  Project  “Im¬ 
proving  Wave  Predictions  Using  Data  Assimilation  in 
a  Spectral  Wave  Model,”  Program  Element  0602435N. 
We  thank  David  Walker  for  graciously  allowing  us  to  use 
his  codes  for  computing  conjugate  gradients  and  tracking 
the  convergence  of  the  cost  function  in  the  adjoint  model, 
and  several  anonymous  reviewers  for  their  many  helpful 
comments  and  suggestions. 


ym=N(Xm’tm)  +  em’  1  (A2) 

in  which  sm  represents  measurement  error.  These  ob¬ 
servations  are  collected  at  discrete  points  distributed  in 
space  and  time.  An  observation  operator  [e.g.,  Hk  in 
Eq.  (3)]  will  not  be  utilized  in  this  derivation  in  accor¬ 
dance  with  the  preceding  study,  for  which  measurements 
are  collocated  with  domain  grid  points  and  the  observed 
quantity  is  the  wave  action  N  rather  than  another  de¬ 
rived  function  of  N.  Because  of  various  sources  of  error, 
the  model  trajectory  generally  deviates  from  the  ob¬ 
servations.  The  role  of  data  assimilation  is  to  minimize 
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the  aggregate  discrepancy  over  the  given  time  interval,  for  which  a  reliable  forecast  will  be  made.  The  aggregate 
It  is  assumed  that  if  the  said  discrepancy  is  minimized,  model-data  discrepancy  is  represented  in  a  generalized 
then  the  optimal  trajectory  will  yield  a  final  state  N(x,  T)  cost  function  of  the  form 


L  cT  rL  rL  rL 

f(x ,  t)Wf(x ,  t ,  x',  t')f(x tf)  dx'  dt'  dx  dt  +  i(x)Wt(x ,  x')i(x')  dxf  dx 

o  Jo  Jo  7  Jo  Jo 

(•r  mm 

+  b(t)Wb{t,t’)b{t’)dt'dt+  X  I  \ym-N(xm,tm)]R^\yn-N(xn,tn)],  (A3) 

JO  JO  m  =  l  n—1 


where  Wj,  Wb  and  Wb  are  the  respective  inverses  of  Cj, 
Q,  and  Q,  defined  in  the  convolution  sense: 

|  |  Wf(x,t,x',t')Cf(x',t',x'',t")  dx'dt'=8(x  —  x")S(t  —  t") 

[  W.(x,x')C.(x' ,x")dx'  =  8(x-x") 

Jo 

[  Wb(t,t')Cb(t',t")dt'  =  8(t-t") 

Jo 

(A4) 

with  the  unit  impulse  8(z)  =  1  for  z  =  0  but  otherwise  zero, 
and  Rmn  are  the  components  of  R,  the  observation  error 
covariance;  that  is,  R  =  eeT,  with  the  overbar  denoting  the 
average  in  the  ensemble  sense.  It  is  important  to  note  that 
the  cost  function  (A3)  is  quadratic  in  the  errors;  thus,  J  has 
one  minimum.  For  the  sake  of  clarity  and  without  loss  of 
generality,  it  is  assumed  that  the  errors  are  uncorrelated 
with  constant  variance.  The  cost  function  becomes 


J  =  W4 


f(x,t)  dxdt  +  Wt 


o 


i(xfdx  +  Wb 


M  X  [ym-N(xm,tm)]2, 

m= 1 


b(t)  dt 


(A5) 


where  co  is  the  inverse  of  the  data  error  variance  [for  the 
steady-state,  homogeneous  form  of  (Al),  Eq.  (A5)  reduces 
to  the  form  of  section  2,  Eq.  (3)].  After  some  standard 
manipulation,  in  which  we  introduce  the  adjoint  variable 


A  =  W/ 


dN  dN 

—  +  c - 

dt  dx 


(A6) 


The  equations  in  (A6)  and  (A7)  constitute  the  Euler- 
Lagrange  conditions  for  local  extrema  of  the  cost 
function. 

While  there  are  several  techniques  for  obtaining  a  so¬ 
lution  to  (A6)  and  (A7),  the  present  study  makes  use  of 
the  conjugate  gradient  method  detailed  in  Walker  (2006) 
and  Polak  and  Ribiere  (1969),  in  combination  with  the 
method  of  representers.  The  optimal  solution  or  “best 
estimate”  N  is  obtained  by  iteratively  solving  the  above 
equations,  reorganized  into  the  coupled  system 

dA  dA  ^ 

- — -c  —  =  -co  y[ N(x^dr^-y™]8(x-xJ\8(t-t) 


(B){ 


A(L,f)  =  0 


{X(x,T)  =  0 


(A8) 


dN  dN  „  TI7_, , 

—  +  c - F  +  Wf  A 

dt  dx  7 

N(x,0)=I+  ' 

[N(0,t)  =  B(t)  +  cW^\(0,t) 


(A9) 


Note  that  we  now  also  have  the  best  estimate  of  the 
errors  /,  /,  and  b  as  the  second  terms  in  all  three  equa¬ 
tions  of  (A9),  respectively.  The  representer  method 
uncouples  the  above  system  by  introducing  representer 
functions  rm(x,t),  1  <  m  <  M.  Each  representer  function 
has  an  adjoint  am(x ,  t)  that  satisfies 


we  obtain 


dA  dX  ^  rx7/ 

x8(x-xm)S(t-tm)  =  0, 

A(L,0  =  0, 

A(x,r)  =  0, 

-c\(0,t)  +  Wb[N(fi,t)~B(t)]  =  0,  and 
-A(jc,0)  +  W.[A(x,0)-/(jc)]  =  0.  (A7) 


OBJ 


da  da 

__m  _c_m=8(x_  x  {  J 

dt  dx  V  m’  V  m> 

%(0()  =  0 

am(x’  r)  =  0 


(A10) 


d!jn+cdIm  =  W- 1, 
dt  dX  * 


(PJ{ 


rJx,0)  =  Wr^am(x,0)- 
Jm(0,t)=cW^aJ0,t) 


(All) 
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Next,  we  seek  the  solution  N  of  the  form 

M 

N(x,  t)  =  Nf(x,  0  +  I  pmrm(x, t) ,  (A12) 

m— 1 

where  the  coefficients  /3m  are  unknown  constants  and 
Np(x,  t )  is  the  prior  estimate  obtained  as  the  solution  of 
(Al)  without  the  errors: 

dN  dN 
—  T"  c —  =  F 
dt  dx 

N(x,0)  =  I(x) 

N(0,t)  =  B(t).  (A13) 

We  combine  (All),  (A12),  and  (A13)  to  get 

M 

D(N)  =  D(Nf)+  I  PmD[rJx,t)\ 

m— 1 
M 

=  F+Wf1  I  Pmam(x,t),  (A14) 

m— 1 


The  optimal  representer  coefficients  still  depend  on  the 
optimal  solution  N ,  so  further  manipulation  is  required. 
Substituting  N  in  (A12)  into  (A17)  gives 


P,n  =  ~0J 


M 


NF(Xm •  O  +  I  /VAm’  O  -  yn 


1=1 


M 


-0)[  NFm  +  &irlm  ~y„ 

1=1 


(A18) 


Thus, 


M 

I 

1=1 


+  W  ym  AFm,  (A19) 


where  5/?m  is  the  Kronecker  delta.  In  matrix  notation, 
the  M  equations  in  (A19)  for  the  representer  coefficients 
Pm  are 


(R  +  fc>  1 1)/3  =  h  =  y  Nf  .  (A20) 


in  which  D  =  (d/dt)  +  c(d/dx ).  Now  by  definition  of  A  and  The  optimal  solution  is  therefore  obtained  as 

by  virtue  of  (A14), 

M  N(x,t)  =  NF(x,t)  +  (y-NF)T(R  +  io~1\y1r(x,t).  (A21) 

A  =  Wf[D(N)-F]  =  I  Pmam(x,t).  (A15) 

m = 1 


Applying  the  differential  operator  D  on  (A15)  and  using 
the  first  equation  of  (A8), 


APPENDIX  B 


M  M 

~D( A)=-  I  PmD(am)=  I  PmS(x-xJS(t-tJ 

m=l  m  =  1 

M 

=  -«  I  [N (xmFj  -  ym]8(x  -  Xm)8(t  -  tj . 

m=l 

(A16) 

Equating  the  coefficients  of  (3m  in  (A16)  yields  the  op¬ 
timal  choice  of  the  representer  coefficients: 

Pm  =  -<o[N(xm,tJ-ym],  1  (A17) 


Discretized  Equations 


In  expanded  form,  the  stationary  homogeneous  wave 
action  equation  [Eq.  (2)]  may  be  written  as 


dC^N 

dx 


+ 


dy 


dCN 

da 


dCnN 
— =  0, 
de 


(Bl) 


where  all  terms  have  been  defined  in  section  2.  Dis¬ 
cretized  using  a  first-order  backward  difference  scheme 
for  spatial  dimensions  and  a  central  difference  for 
spectral  dimensions,  this  equation  becomes 


1  1 

\  r  (  CxKjJXl  ^xi- 1  ~  1  J’k,l  )  Ay(CyjNW  Cyj_]  ^i,j—l,k,l^ 

+ i[(1 "  + -  d  ■ +  JV- 1,/] 

+  We[{l~  ’Wi  +  -  (1  +  v)CeilNim_  i]  =  0, 


(B2) 


where  indices  i, A,  and  /  track  the  spatiospectral  grid  are  expressed  as  weighted  central  differences  with 

location  in  terms  of  coordinates  jc,  y,  cr,  and  6 ,  respec-  weights  ^  and  77  between  zero  and  one.  Taking  the 

tively.  Here,  the  frequency  and  directional  derivatives  adjoint  of  (B2)  gives 
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Cxp^i,j,k,l  ^ i+lj,k,l )  +  Ay  Cyp'-ilKl  ^ij+ljcj)  +  2 Act CT'  V^i,j,k-\,l  +  ^v^i,j,k,l  0  +  V^ij,k+A 

+  WecA  -  1 +  2^im  ~  (! +  ^m+i\  =  (ym  ~  N Jsii  Asisi  • 


(B3) 


where  Xqxi  *s  defined  as  the  residual  variable  or  in¬ 
novation  (the  difference  between  model  estimate  and 
observation)  at  the  spatiospectral  location  given  by  in¬ 
dices  i,  j,  k ,  and  /,  and  the  observation  locations  are  de¬ 
noted  by  m,  where  1  <  m  <  M.  Each  8  is  a  Kroenecker 
delta  for  a  given  spatiospectral  index,  which  is  equal  to 
one  only  at  the  index  of  an  observation  location  and  is 
zero  otherwise.  Note  that  in  (B3),  the  energy  propa¬ 
gation  speeds  Cx,y^6  are  assumed  to  be  locally  con¬ 
stant  with  respect  to  coordinates  x,  y,  cr,  and  0 ,  which 
leads  to  a  different  form  for  the  coefficients  in  (B3) 
compared  to  those  in  (B2).  The  adjoint  to  linear  ho¬ 
mogeneous  SWAN  essentially  solves  (B3),  forced  at 
selected  interior  grid  observation  locations  [x(/m), 
y(jm)]  with  innovations  (ym  -  Nm),  each  of  which  in¬ 
cludes  observed-minus-estimated  spectral  density  for 
frequencies  cr(km)  and  directions 
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